#include <cmath>

#include "sunposition.h"


// 计算正弦函数，结果保留指定精度
double SunPosition::MathSin(double sin, int precision)
{
    sin = round(sin * 1e5) / 1e5; // 保留5位小数
    return round(std::sin(sin) * pow(10, precision)) / pow(10, precision);
}

// 计算余弦函数，结果保留指定精度
double SunPosition::MathCos(double cos, int precision)
{
    cos = round(cos * 1e5) / 1e5; // 保留5位小数
    return round(std::cos(cos) * pow(10, precision)) / pow(10, precision);
}

// 计算时区
int SunPosition::calculateTimeZone(double longitude)
{
    int cutValue = static_cast<int>(longitude / 15);
    double modeValue = std::abs(longitude) - cutValue * 15;
    if (modeValue <= 7.5)
        return cutValue;
    else
        return cutValue + (longitude > 0 ? 1 : -1);
}

// 计算太阳位置，返回高度角和方位角
std::pair<double, double> SunPosition::calculateSunPosition(double longitude, double latitude, int year, int month, int day, int hour, int minute, int second)
{
    double A = year / 4.0;
    double B = A - std::floor(A);
    double C = 32.8;
    if (month <= 2)
        C = 30.6;
    if (B == 0.0 && month > 2)
        C = 31.8;
    double G = std::floor(30.6 * month - C + 0.5) + day;
    double L = longitude / 15.0; // 经度修正

    int timeZone = calculateTimeZone(longitude);

    double H = hour - timeZone + minute / 60.0; // 时刻修正
    double N = G + (H - L) / 24.0; // 计算积日
    double theta, NO; // theta为日角
    NO = 79.6764 + 0.2422 * (year - 1985) - std::floor((year - 1985) / 4.0);
    theta = 2 * M_PI * (N - NO) / 365.2422;

    double lat; // 赤纬角
    lat = 0.3723 + 23.2567 * MathSin(theta) + 0.1149 * MathSin(2 * theta) - 0.1712 * MathSin(3 * theta) - 0.758 * MathCos(theta) + 0.3656 * MathCos(2 * theta) + 0.0201 * MathCos(3 * theta);
    lat = lat * M_PI / 180.0;

    double Et; // 时差（以分钟为单位）
    Et = 0.0028 - 1.9857 * MathSin(theta) + 9.9059 * MathSin(2 * theta) - 7.0924 * MathCos(theta) - 0.6882 * MathCos(2 * theta);

    double timeAngle; // 时角
    timeAngle = ((hour - 12) + (minute - (120.0 - longitude) * 4.0 + Et) / 60.0) * 15.0;
    timeAngle = timeAngle * M_PI / 180.0;

    // 将纬度由角度转换为弧度
    latitude = latitude * M_PI / 180.0;

    // 计算太阳高度角
    double height = MathSin(latitude) * MathSin(lat) + MathCos(latitude) * MathCos(lat) * MathCos(timeAngle);
    height = std::asin(height);

    // 计算太阳方位角
    double direction = (MathSin(height) * MathSin(latitude) - MathSin(lat)) / (MathCos(height) * MathCos(latitude));
    direction = std::acos(direction);

    if (timeAngle < 0) // 以正南方向为0度，顺时针为正
        direction *= -1.0;

    return { height, direction };
}
